Driving defect modes of Bose-Einstein condensates in optical lattices 
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We present an approximate analytical theory and direct numerical computation of defect modes 
of a Bose-Einstein condensate loaded in an optical lattice and subject to an additional localized 
(defect) potential. Some of the modes are found to be remarkably stable and can be driven along 
the lattice by means of a defect moving following a step-like function defined by the period of 
Josephson oscillations and the macroscopic stability of the atoms. 
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Introduction.- Loading and manipulating cold bosonic 
atoms in optical lattices is a fascinating, rapidly growing 
branch of cold atom physics [B 0- Q Recent progresses 
of experimental investigations of Bose-Einstein conden- 
sates (BECs) in optical lattices (OLs) |3( and in par- 
ticular the direct observation of a gap soliton p| have 
stimulated studies of the nonlinear dynamics of matter 
waves in periodic media (see e.g. [2l LJ and references 
therein). OLs provide a tool for changing the effective 
properties of atomic media allowing for the existence of 
solitary waves in a condensate with positive scattering 
length. Moreover, OLs provide also ways for manipulat- 
ing nonlinear matter waves. For instance, smooth mod- 
ulations of an OL can be used to control the dynamics 
of gap solitons On the other hand, strongly local- 

ized defects generated by narrow laser beams have been 
used in experiments to generate rotational motion of con- 
densed atoms |j| and to study scattering of a soliton on 
a defect for the sake of understanding BEC expansion in 
an effectively one-dimensional (ID) random potential ||. 

In this Letter we study the combined effect of an OL 
and of a defect, whose spatial extension is of the order of 
the OL period, on the spatial distribution and dynamics 
of matter waves. Specifically we study the existence and 
stability of stationary modes localized in the vicinity of 
the defect, concentrating on the approximate analytical 
description of defect modes and the numerical simulation 
of their dynamics, establishing the stability properties of 
the modes. We also show that a defect mode can be 
driven through an OL over hundreds of lattice periods, 
thus representing an effective tool for matter waves man- 
agement. 



(rj/y/2n£) exp (— a; 2 /2£ 2 ), is a defect potential character- 
ized by its amplitude 77, width £ < tt (i.e. localized on 
a distance of the order of one lattice period) and center 
Xd- This potential can be experimentally obtained by 
illuminating the condensate with an orthogonal narrow 
Gaussian laser beam. 

Since the defect is strongly localized we can use the 
expansion over the orthonormal basis of Wannier func- 
tions (WFs), which are also localized on about one lat- 
tice period |9j- To do this we first expand the solution 
of JU, tp(x,t) = J2a=oJ dkif ak (x)f a (k,t), over the set 
of Bloch functions (BFs) (f a k{x) solving the eigenvalue 
problem (— d 2 /dx 2 + Acos(2x)) ip a k — E a (k)ip a k where 
a = 0, 1, ... stands for the number of the zone, k is the 
wave vector considered in the first Brillouin zone (BZ), 
k E [— 1 , f ] , and hereafter integrals with respect to k are 
over the first BZ. Substituting the above expansion in 
Eq. 0) one obtains that the envelope f a (k) satisfies 



if a (k)-E a (k)f a (k) 



dxip ak (x)\ljj(x, t)flp(x, t) 



dk' / dxtp ak (x)V d (x - x d )f a 'k'{x)f a '{k'). (2) 



The model.- In the quasi ID limit a BEC is governed b 
the dimensionless ID Gross-Pitaevskii (GP) equation 



iipt = -ip xx + A cos(2x)-0 + V d (x~ x d )ip + o-\ip\ 2 ip, (1) 

where a = 1 for repulsive and a = — 1 for attrac- 
tive interactions among atoms, A is the amplitude of 
the lattice potential, the lattice constant is chosen to 
be 7T without restriction of generality and Vd(x) — 



The integration with respect to x is carried out over the 
real axis hereafter and to shorten notations we drop the 
temporal argument of the envelope f a . 

Next we use the standard definition of WFs: w an (x) = 
(1/V2) / fak(x) exp(—iTrnk)dk which implies <p a k(x) — 
(1/V2) w an{x) exp(i7rn/c) and introduce the matrix 
element V™ 1 = (1/2) / dxVd(x)w an (x)w a ' n '(x) describ- 
ing probability of transitions between two bands, a and 
a' and between sites n and n' as well as its Fourier 
transform V aoi >(k,k') = Y, n ,n> V aa< ex P {^(k'n 1 - kn)). 
Then, the last term on the right hand side of J2J is rewrit- 
ten in the form Y^ Q , J dk'V aa '(k, k')f a {k') . 

Approximate description.- Let us start with the case 
where the defect is localized at the origin, Xd = 0, which 
also coincides with a minimum of the lattice potential, 
A < 0. Matter-wave modes of small amplitude, having 
width A, exceeding a lattice period, i.e. satisfying A 3> vr, 



2 



are well described in the effective mass approximation 
(EMA) 0. In that case, the characteristic scale of f a (k) 
in the momentum space is much smaller than the vector 
of the reciprocal lattice which in our case is equal to 2, 
i.e. Afc <C 1, the latter being the scale of variation of 
V aa '(k,k'). In other words, in the EMA the envelope 
f a (k) is much more localized than V aa i(k, k'). 

Let us consider localized modes whose carrier wave 
wavevector fco borders the BZ: ko = 1. Then the spa- 
tial domain where f a (k) significantly differs from zero is 
given by k G (1 — Afc, 1 + Afc) allowing us to approximate 
with exponential accuracy 



V aa ,(k,k')f a ,(k')dk' 



*aa f 

~2~r 



f a .{k')dk'. (3) 



Hereafter V aa ' = 27rV r Q , ce / (1, 1), V a = V aa . 

In order to simplify the nonlinear term in we use 
the properties of the WFs (see e.g. [ljj) an d rewrite it 
as follows 



\ E E / dkidk 2 dk 3 S k+klM+k3 W°ZZZ 



rOnin 2 n 3 
*3 

^1,2,3 ^1,2,3 ' 

piw (fe 2 „ 2 +fe3„ 3 -/ S1 n 1 ) 7ai {kl)fa2 {k2)fa3 { f, 3) (4) 



xe 



where W™l2lal = I dxw an w aini w a2n2 w a3n3 . 

The EMA is justified when the effective mass, M a = 
[d 2 E a (k)/dk 2 ]~ , is of order one and the respective WFs 
are localized on a very few lattice periods. These condi- 
tions are verified for OLs with amplitudes of order of a 
few recoil energies (in dimensionless units for \A\ ~ 1) 
what is illustrated by examples in Table [I] 

Subject to the mentioned conditions we rewrite (@J) as 

Zjai 2 3 Waa 1 a 2 a 3 Qa 1 a 2 a 3 where 

M7 — ( T\n2+n 3 -n 1 Tjr0n 1 n 2 n 3 (r\ 

VVaa ia2 a 3 — ~Z 2_^i V > VV a ai a 2 a 3 ! \° ) 



dkidk 2 f a (ki + k 2 - fc)/ ai {k\)f a2 (k 2 ). 



Let us now take into account that for our lattice poten- 
tial the WFs have the symmetry w an (x) = w an (—x) for 
a = 0, 2, ... and w an (x) = —w an (—x) for odd a. Then 
one readily concludes that the even defect Vd (x) results in 
inter-band transitions only among the bands having the 
same parity what is reflected by the property V aa i = 
whenever \a — a'\ is odd. On the other hand the nonlin- 
earity, while couples different modes, due to Woi, does 
not affect the particle exchange among them preserving 
the averaged densities. Thus, assuming that initially all 
atoms were loaded in one of the lowest bands one can em- 
ploy the tight-binding approximation which neglects the 
inter-band transitions and results in decoupling of the 
equations for the envelopes f a (k). Then Eq. J2J become 

ifa{k) = E a {k)f a {k) + V a f a {k)dk + oW aa Q aaa . 



Near the boundary of the BZ we can expand E a (k) w 
E a + (k — k ) 2 /2M a where E a — E a (l) is the energy of 
the a-th band at the boundary of the BZ (see Fig. QJ. 
Then, introducing f a (x,t) = e" 4 ^* / dke< k - ko > f a (k,t) 
and after straightforward algebra we obtain the equation 
the nonlinear Shrodinger (NLS) equation with a delta 
impurity 



.df a 

'~dt 



1 d 2 f Q 
2M a dx 2 



<rW aa \f a \ 2 f a = V a 5(x)f a . (6) 



TABLE I: Parameters of the spectrum at the boundary of 
the BZ for the two lowest bands. The integrals W aa ' ar e 
approximated with accuracy of about 3% for A = — 1 by 
considering only neighbouring sites in (JKJ. The accuracy of 
this approximation grows with the OL amplitude. 
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M 


E 


W 00 


Mi 




Woi 


-1 

-5 


-0.163 
-2.5 


0.471 
-2.076 


0.25 
0.358 


0.1 
0.3 


1.467 0.174 
2.5 0.234 


0.182 
0.274 



In Table [I] we present the numerical values of the coef- 
ficients of 0, used below for the numerical simulations. 
The wave function ip(x,t) is recovered from Eq. jnj) by 



^E 



(-l) n w an (x)f a (nL). 



(7) 



Defect modes. - Stationary defect modes of ©have the 



form: f a (x) 



4> a {x) with <f> a (x) real 



® ha 



The 



applicability of the EMA and the requirement for the 
frequency of the mode to belong to a forbidden gap, imply 
the conditions: |e a | <C E a+ i—E a and e a M a < (Fig-QJ. 

We discuss below defect modes in the first lowest gap 
(i.e. for a = 0, 1). For aM a < one obtains cosh-modes 



<t>a{x) 



^2\e a \/W aa 



cosh( v /2|M Q e Q | (\x\ - x a ) 



(8) 



Here x a 
and = 



atanh 



(sign(M a V a )y/e*/e a } /^2\M a e a \, 

-M a V 2 /2. This mode exists when |e Q | > 
For M a V a < the cosh- mode has only one maximum 
and otherwise it has a two- hump profile |ll| . 

When crM a > there is another localized solution of 
Eqs. 10 - a sinh-mode - which corresponds to a smaller 
detuning |e a | < \e*\ such that e a a < 



(f> a (x) 



^/2\e a \/W a 



sinh 



2\M a e a \(\x\ 



(9) 



Now x a — atanh {^/e a /e^j / y/2\M a e a \. 

Since w\ n are odd functions and the defect potential 
is localized on its zero we have that |Vi| *C |Vo| {Vi = 
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M, > 



two-hump cosh-mode 
(a<0) 




cosh-modc (a > 0) — £i 
sinh-mode (a < 0) 



M < 




sinh-mode (a > 0) 

8*| 

£, -\- cosh-mode (a < 0) 

two-hump cosh-mode 
(a>0) 



K,>0 



K<0 



FIG. 1: Location of the impurity modes for Vd > 0, Vd < 0. 



in the limit t — 0). Thus /o is subject to much larger 
influence of the defect than f% . We also notice that since 
M a =0(1) one has |e»| ~ V„. On the other hand the 
applicability of our approach implies that all the terms 
in Eq. © are of the same order and that A ^> tt. Thus 
\e a \ ~ \V a \/X ~ 1/A 2 < \E a \ ~ 1 and V a < 1. Below we 
concentrate on the cosh-mode excited in the vicinity of 
the lowest band (a = 0) and with a = 1. 
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FIG. 2: Approximate analytical (thick lines) and numerical 
(thin lines) shapes of the stationary cosh-mode of the lowest 
band (a,b) A = -1, t = 0.1; (d,e) A = -5, I = 0.1; (c) and 
(f) dynamics of the modes (a) and (d), correspondingly. 

In Fig. |3 we plot both our analytical approximation 
for cosh-modes given by Eqs. J7J) and JSJ) and the direct 
numerical solutions of Eq. (JTJ for different OL and defect 
parameters showing the good accuracy of our one-band 
approximation. We have checked that decreasing the de- 
fect width by a factor of ten, shapes of the defect modes 
do not change and their amplitudes differ only by 1%, 
what confirms the excellent accuracy obtained from the 
delta-approximation for t -C 7r. The analytical approxi- 
mation becomes worse for increasing values of the defect 
strength. This occurs because such a defect originates 
modes localized on a very few lattice periods, and the 
tight-binding approximation should be used instead |ldl | . 

The first simple test of the stability is shown in Fig. [21 
where we have computed the evolution of the approxi- 
mate analytical cosh-modes, which introduces a pertur- 
bation of the order of a few percents with respect to the 



true mode. Next we impose stronger perturbations, by 
shifting the position of the defect. First, we shift the 
defect at time t = 10 to the nearest minimum of the pe- 
riodic potential, = tt, and the mode follows the defect 
displacement. After some relaxation time the profile of 
the defect mode centered about x — tt is very similar to 
the initial one centered about x — [Figs.|!3Vb)]. When 
the defect is shifted to Xd = tt/2, the mode is destroyed, 
since the new position does not correspond to a stable 
configuration [Figs. 02c- d)]. 




FIG. 3: [Color online] (a-b) Density plots of the dynamics 
of the cosh-mode with initial parameters as in Fig. |5Ja) and 
Xd = 0. At t — 10 the defect is shifted to the position (a) 
Xd — tt and (c) Xd = tt/2 indicated by the arrows. In (b) and 
(d) dashed and solid lines correspond to the initial (t = 0) and 
final (t — 100) profiles of the density of the impurity mode. 

A driven defect mode.- The robust mode behavior 
shown in Fig. OJa) suggests that it could be driven by 
moving the defect along the lattice. On the other hand 
the mode is destroyed when it stays close to the local 
maxima of the potential. This suggests that the defect 
motion should be defined as a step-like function com- 
posed of two characteristic time intervals: a fast one r, 
in which the defect is shifted by one lattice period and 
a slow one T, allowing for the cosh-mode to recover its 
shape on the new site. 

In Fig.[S](b) we see that the mode dynamics is reflected 
in the change of its amplitude, what can be interpreted 
as a tunneling of a part of atoms in a time To, attracted 
by the shifted defect. Since this is a "single" tunnel- 
ing process and the intermediate states are unstable the 
characteristic times should satisfy r r <T. 

To estimate the Josephson tunneling time tq we take 
into account that the EMA implies relatively small differ- 
ences in populations of neighbor potential wells and that 
we are near the boundary of the BZ, what means that 
the phase difference between atoms in adjacent cells is n, 
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thus we can use the results of [12J for the estimation of 
the half-period of the oscillations of the linear tunneling 
of a BEC in a double-well potential. This gives r w 3.67 
for A = — 1. Hence, the defect must be shifted faster 
than tt/to ps 0.86 and then stop for T ^> 3.67 to let the 
system relax to equilibrium. 

In Fig. 0] we show numerical simulations of the driv- 
ing of the numerically found defect modes with Xd(t) = 
TrJ2j=o@(t ~ *o — 100 j), where 9(t) is the Heaviside 
step function and to is a time at which motion is 
started. Shown are also the mean dispersion or the 
width of the wavepacket D = \J (x 2 ) — (x) 2 where (x n ) — 
Jx n \if>\ 2 dx/ f\i>\ 2 dx [Fig.Etb-c)]. 




FIG. 5: Adiabatic excitation of the cosh-mode. Shown are 
the initial envelope soliton (dashed line) for A — —1, so = 0.01 
and n = and the final cosh-mode (solid line) corresponding 
to I = 0.1, n = 0.1, with e* = 0.0034 and e Q ~ 0.027. 




FIG. 4: [Color online] (a) Pseudocolor and amplitude plots of 
the motion of the driven defect mode. Initial parameters as 
in Fig. 13 (b) Defect position, Xd(t) (thin line), and average 
coordinate of the center of the defect mode (x) (t) (thick line) . 
(c) Evolution of the dispersion D(t). 

In Fig. 0] (c) we see that the dispersion of the wave 
packet grows in the interval t £ (200, 600) after which 
the moving mode achieves a stationary profile. This re- 
shaping occurs through the emission of radiation which 
also contributes to the appearance of small oscillations of 
(x) [see Fig. 0] (b)] in this time interval. 

Excitation of defect modes.- A way to excite defect 
modes in real situations is to create a gap soliton with 
energy near the band edge in a homogeneous OL , then 
by increasing adiabatically the intensity of the transverse 
laser bean one could generate the defect mode. We have 
realized this numerically (see Fig. [SJ starting from an en- 
velope soliton with energy Eq +e and increasing rj. The 
so obtained final state is the cosh-mode. 

Conclusion.- We have shown that a BEC in a ID OL 
with an additional localized potential supports defect 
modes. In the weak density limit such modes are well 
described by the effective mass approximation reducing 



the problem to the NLS equation with a delta impurity. 
We have constructed the defect modes analytically and 
numerically and have shown that some of such modes 
can be driven along the lattice by a defect which follows 
a specific law of motion. Finally we have shown that de- 
fect modes can be easily excited by adiabatically switch- 
ing on the defect. Our results open new possibilities for 
controlling matter waves which we hope will stimulate 
further experiments with Bose-Einstein condensates in 
optical lattices. 
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